home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / rng / knuthran.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-12-02  |  4.5 KB  |  178 lines
  1. /* rng/knuthran.c
  2.  * 
  3.  * Copyright (C) 2001 Brian Gough, Carlo Perassi
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /*
  21.  * This generator is taken from
  22.  *
  23.  * Donald E. Knuth
  24.  * The Art of Computer Programming
  25.  * Volume 2
  26.  * Third Edition
  27.  * Addison-Wesley
  28.  * Section 3.6
  29.  *
  30.  * The comments are taken from the book
  31.  * Our comments are signed
  32.  */
  33.  
  34. #include <config.h>
  35. #include <stdlib.h>
  36. #include <gsl/gsl_rng.h>
  37.  
  38. #define BUFLEN 2009        /* [Brian]: length of the buffer aa[] */
  39. #define KK 100            /* the long lag */
  40. #define LL 37            /* the short lag */
  41. #define MM (1L << 30)        /* the modulus */
  42. #define TT 70            /* guaranteed separation between streams */
  43.  
  44. #define evenize(x) ((x) & (MM - 2))    /* make x even */
  45. #define is_odd(x) ((x) & 1)    /* the units bit of x */
  46. #define mod_diff(x, y) (((x) - (y)) & (MM - 1))    /* (x - y) mod MM */
  47.  
  48. static inline void ran_array (unsigned long int aa[], unsigned int n,
  49.                   unsigned long int ran_x[]);
  50. static inline unsigned long int ran_get (void *vstate);
  51. static double ran_get_double (void *vstate);
  52. static void ran_set (void *state, unsigned long int s);
  53.  
  54. typedef struct
  55. {
  56.   unsigned int i;
  57.   unsigned long int aa[BUFLEN];    /* [Carlo]: I can't pass n to ran_array like
  58.                    Knuth does */
  59.   unsigned long int ran_x[KK];    /* the generator state */
  60. }
  61. ran_state_t;
  62.  
  63. static inline void
  64. ran_array (unsigned long int aa[], unsigned int n, unsigned long int ran_x[])
  65. {
  66.   unsigned int i;
  67.   unsigned int j;
  68.  
  69.   for (j = 0; j < KK; j++)
  70.     aa[j] = ran_x[j];
  71.  
  72.   for (; j < n; j++)
  73.     aa[j] = mod_diff (aa[j - KK], aa[j - LL]);
  74.  
  75.   for (i = 0; i < LL; i++, j++)
  76.     ran_x[i] = mod_diff (aa[j - KK], aa[j - LL]);
  77.  
  78.   for (; i < KK; i++, j++)
  79.     ran_x[i] = mod_diff (aa[j - KK], ran_x[i - LL]);
  80. }
  81.  
  82. static inline unsigned long int
  83. ran_get (void *vstate)
  84. {
  85.   ran_state_t *state = (ran_state_t *) vstate;
  86.  
  87.   unsigned int i = state->i;
  88.  
  89.   if (i == 0)
  90.     {
  91.       /* fill buffer with new random numbers */
  92.       ran_array (state->aa, BUFLEN, state->ran_x);
  93.     }
  94.  
  95.   state->i = (i + 1) % BUFLEN;
  96.  
  97.   return state->aa[i];
  98. }
  99.  
  100. static double
  101. ran_get_double (void *vstate)
  102. {
  103.   ran_state_t *state = (ran_state_t *) vstate;
  104.  
  105.   return ran_get (state) / 1073741824.0;    /* [Carlo]: RAND_MAX + 1 */
  106. }
  107.  
  108. static void
  109. ran_set (void *vstate, unsigned long int s)
  110. {
  111.   ran_state_t *state = (ran_state_t *) vstate;
  112.  
  113.   long x[KK + KK - 1];        /* the preparation buffer */
  114.  
  115.   register int j;
  116.   register int t;
  117.   register long ss = evenize (s + 2);
  118.  
  119.   for (j = 0; j < KK; j++)
  120.     {
  121.       x[j] = ss;        /* bootstrap the buffer */
  122.       ss <<= 1;
  123.       if (ss >= MM)        /* cyclic shift 29 bits */
  124.     ss -= MM - 2;
  125.     }
  126.   for (; j < KK + KK - 1; j++)
  127.     x[j] = 0;
  128.   x[1]++;            /* make x[1] (and only x[1]) odd */
  129.   ss = s & (MM - 1);
  130.   t = TT - 1;
  131.   while (t)
  132.     {
  133.       for (j = KK - 1; j > 0; j--)    /* square */
  134.     x[j + j] = x[j];
  135.       for (j = KK + KK - 2; j > KK - LL; j -= 2)
  136.     x[KK + KK - 1 - j] = evenize (x[j]);
  137.       for (j = KK + KK - 2; j >= KK; j--)
  138.     if (is_odd (x[j]))
  139.       {
  140.         x[j - (KK - LL)] = mod_diff (x[j - (KK - LL)], x[j]);
  141.         x[j - KK] = mod_diff (x[j - KK], x[j]);
  142.       }
  143.       if (is_odd (ss))
  144.     {            /* multiply by "z" */
  145.       for (j = KK; j > 0; j--)
  146.         x[j] = x[j - 1];
  147.       x[0] = x[KK];        /* shift the buffer cyclically */
  148.       if (is_odd (x[KK]))
  149.         x[LL] = mod_diff (x[LL], x[KK]);
  150.     }
  151.       if (ss)
  152.     ss >>= 1;
  153.       else
  154.     t--;
  155.     }
  156.  
  157.   state->i = 0;
  158.  
  159.   for (j = 0; j < LL; j++)
  160.     state->ran_x[j + KK - LL] = x[j];
  161.   for (; j < KK; j++)
  162.     state->ran_x[j - LL] = x[j];
  163.  
  164.   return;
  165. }
  166.  
  167. static const gsl_rng_type ran_type = {
  168.   "knuthran",            /* name */
  169.   0x3fffffffUL,            /* RAND_MAX *//* [Carlo]: (2 ^ 30) - 1 */
  170.   0,                /* RAND_MIN */
  171.   sizeof (ran_state_t),
  172.   &ran_set,
  173.   &ran_get,
  174.   &ran_get_double
  175. };
  176.  
  177. const gsl_rng_type *gsl_rng_knuthran = &ran_type;
  178.